/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Functions to perform generic operations on the image arrays.

// $Id$

#include "imops.h"
#include "blitz-fits.h"
#include "fits-utilities.h"
#include "vecops.h"
#include "mcrx-debug.h"

using namespace blitz;
using namespace CCfits;

mcrx::array_1 
mcrx::integrate_image(CCfits::FITS& file, 
		      const std::string& hdu_name,
		      T_float distance_factor,
		      const std::string& column_name,
		      const array_1& wavelengths,
		      const std::string& bol_keyword,
		      const T_unit_map& units,
		      boost::function<terminator::t_type()> term)
{
    Array< float, 3> image;

    try {
      ExtHDU& hdu = open_HDU (file, hdu_name);
    }
    catch (CCfits::FITS::NoSuchHDU&) {
      cout << hdu_name << " not found, skipping" << endl;
      array_1 dummy(wavelengths.size());
      dummy=0;
      return dummy;
    }

    ExtHDU& hdu = open_HDU (file, hdu_name);
    cout << "\nReading " << hdu_name << endl;
    read (hdu, image);
    ASSERT_ALL (image < blitz::huge (T_float() ));

    // If compression was used on the data cubes, the image values can
    // be slightly negative, so an assertion here would strictly not
    // be so good.  Instead we print a warning message
    if (any ( image < 0))
      cout << "  Warning: negative numbers found in image!  Minimum value: " 
	   << min (image) << endl;
    if (term ()) return array_1 ();

    T_float sb_factor;
    hdu.readKey("sb_factr", sb_factor);

    array_1 sed(image.extent(thirdDim));

    cout << "  Integrating image" << endl;

    // Now contract the xy dimensions. In the process we also
    // convert surface brightness to inferred luminosity using the
    // pixel solid angle and camera distance.
    array_2 temp(image.extent(firstDim),
		 image.extent(thirdDim));
    temp =sum (image (tensor::i, tensor::k, tensor::j), tensor::k);
    sed = sum ( temp  (tensor::j, tensor::i), tensor::j)*
      (distance_factor/sb_factor);
    ASSERT_ALL (sed < blitz::huge (T_float ()));
    if (term ()) return array_1 ();

    // Calculate bolometric luminosity
    const T_float L_bol = integrate_quantity (sed, wavelengths);
      
    cout << "  Bolometric luminosity: " << L_bol << endl;
      	
    // write columns to INTEGRATED_QUANTITIES HDU
    ExtHDU& iq_HDU = open_HDU (file, "INTEGRATED_QUANTITIES"); 
    iq_HDU.addKey(bol_keyword, L_bol,
		  "[" + units.get("luminosity") +
		  "] Bolometric luminosity in  " + hdu_name);
    if(column_name != "") {
      try {
	iq_HDU.column(column_name);
      }
      catch (Table::NoSuchColumn&) {
	iq_HDU.addColumn(Tdouble, column_name, 1,
			 units.get("L_lambda"));
      }
      write (iq_HDU.column(column_name), sed, 1);
    }

    return sed;
}
